Parameter Dependent Problems

56. Parameter Dependent Problems#

We have observed that adding the constraint to the A-block may improve the convergence of saddle-point solvers. If we add the constraint with a large factor $1/\varepsilon$, we may consider the augmented problem as penalty formulation for the constrained minimization problem, and skip the second equation and the Lagrange parameter completely. We call this a parameter dependent problem:

\[ \big(A + \frac{1}{\varepsilon} B^T C^{-1} B \big) u = f \]

Let $V_0 = \operatorname{ker} B$ the null-space of the $B$-matrix.

On the null-space $V_0$ only the $A$-operator is seen. On its complement $V_0^\bot$, the second part $\frac{1}{\varepsilon} B^T C^{-1} B$ dominates for small $\varepsilon$. Unless $V_0$ is trivial, the condition number degenerates with $O(\varepsilon^{-1})$.

The goal is to develop preconditioners for the parameter dependent problems which are robust in $\varepsilon$.

56.1. Dirichlet boundary conditions by penalty:#

Adding the Dirichlet condition $u = u_D$ on $\Gamma_D$:

\[ \int \nabla u \nabla v + \int_{\Gamma_D} \alpha u v = \int f v + \int_{\Gamma_D} \alpha u_D v \]

with large parameter $\alpha = \varepsilon^{-1}$, or depending on the mesh-size $\alpha = h^{-1}$.

The null-space $V_0$ is the finite element sub-space of $H_{0,D}^1 = \{ v \in H^1 : v = 0 \text{ on } \Gamma_D \}$.

56.2. Penalty formulation for the flux:#

Approximating the mixed formulation $\sigma = \nabla u$, and $\operatorname{div} \sigma = -f$ by penalty leads to:

\[ \int \sigma \tau + \frac{1}{\varepsilon} \int \operatorname{div} \sigma \operatorname{div} \tau = \frac{-1}{\varepsilon} \int f \operatorname{div} \tau \]

This is a variational formulation on the Hilbert-space $H(\operatorname{div}) := \{ \tau \in [L_2]^d : \operatorname{div} \tau \in L_2\}$.

The null-space $V_0$ consists of all $\operatorname{div}$-free functions.

56.3. Maxwell equations:#

Very similar is the equation for the magnetic field, a special case of the Maxwell equations. Given a divergence-free current density $j$, search for a divergence-free magnetic field $B$ such that

\[ \operatorname{curl} \mu^{-1} B = j, \]

where $\mu$ is the so called permeability, a kind of conductivity for the magnetic field. Since we search for a $\operatorname{div}$-free $B$, and assuming a simply connected domain, we may introduce a vector-potential $u$ such that $B = \operatorname{curl} u$. Together we have the second order equation

\[ \operatorname{curl} \mu^{-1} \operatorname{curl} u = j. \]

Its weak formulation (skipping discussion of boundary conditions) is

\[ \int \mu^{-1} \operatorname{curl} u \, \operatorname{curl} v = \int j v \qquad \forall \, v \]

There is no unique solution: adding an arbitrary gradient field to a solution is also a solution. One way out is Coloumb - gauging, where we add the constraint $\int u \nabla \psi = 0$, i.e. a weak formulation of $\operatorname{div} u = 0$. Another, very simple method is to regularize with a small $L_2$-term:

\[ \int \mu^{-1} \operatorname{curl} u \, \operatorname{curl} v + \int \varepsilon u v = \int j v \]

Again, we have a parameter dependent problem, with the fictitious regularization parameter $\varepsilon$. The canonical space is

\[ H(\operatorname{curl}) = \{ u \in [L_2]^3 : \operatorname{curl} u \in [L_2]^3 \}, \]

equipped with the norm $$ \| u \|_{H(\operatorname{curl})}^2 = \| u \|_{L_2}^2 + \| \operatorname{curl} u \|_{L_2}^2 $$

The null-space is $$ V_0 = \{ v : \operatorname{curl} v = 0 \} = \nabla H^1 $$

56.4. Penalty formulation for the Stokes equation:#

We approximate the div-free condition by penalty:

Find $u \in [H_{0,D}]^d$ such that

\[ \int \nu \nabla u : \nabla v + \frac{1}{\varepsilon} \int \operatorname{div} u \operatorname{div} v = \int f v \]

It looks similar as the penalty formulation for the flux, but now the A-form is a $H^1$-inner product instead of an $L_2$-inner product. Now, conforming low order methods lead to bad approximations (called locking effect). Robust discretizations are constructed and analyzed by by mixed formulations. If the pressure space consists of discontinuous finite elements, and the $C$ matrix is chosen as $L_2$-inner product, it can be cheaply inverted. We end up with the discretization of the parameter dependent problem. However, the detour over the mixed mixed formulation leads to a projection operator $R_h$:

Find $u_h \in V_h \subset [H_{0,D}]^d$ such that

\[ \int \nu \nabla u_h : \nabla v_h + \frac{1}{\varepsilon} \int R_h \operatorname{div} u_h \operatorname{div} v_h = \int f v_h \]

A stable example is to choose $V_h$ a $P^2$ finite element space. Then $\operatorname{div} u_h$ is a $P^1$ function. Choosing a $P^0$ pressure leads to a projection operator $R_h$ which is the $L_2$-orthogonal projection from $P^1$ onto $P^0$.

57. Robust two-level methods for parameter dependent problems#

Let

\[ A_h^\varepsilon = A + \frac{1}{\varepsilon} B^T C^{-1} B \]

be the discretization matrix of the parameter dependent problem on the finite element space $V_h$. We are going to analyze two-level preconditioners

\[ C_{2L}^{-1} = D_h^{-1} + P \, ( A_H^\varepsilon ) ^{-1} P^T, \]

where $A_H^\varepsilon$ is the discretization matrix on the coarse space $V_H$, $P : V_H \rightarrow V_h$ is the prolongation matrix, and $D_h$ is a local block-Jacobi preconditioner on $V_h$.

57.1. Robust smoothers#

By the ASM - lemma, we can represent the norm generated by the smoother $D_h$ as

\[ \| u \|_{D_h}^2 = \inf_{u = \sum u_i} \sum \| u_i \|_{A + \frac{1}{\varepsilon} B^T C^{-1} B}^2, \]

where we have the sub-space splitting

\[ V_h = \sum V_i. \]

If $u$ is a function in the finite-element null-space $V_0$, then its norm $\| u \|_{A + \frac{1}{\varepsilon} B^T C^{-1} B}^2$ is independent of the small parameter $\varepsilon$. Only if we can split $u$ into local null-space functions $u_i$, the $D_h$ norm is also independent of $\varepsilon$. Thus, we want a splitting compatible with the null-space:

\[ V_0 = \sum V_i \cap V_0 \]

Example: Maxwell equations

For a compatible discretization with Nedelec finite elements, the discrete null-space consists of discrete gradients:

\[ \{ v_h \in V_h : \operatorname{curl} v_h = 0 \} = \{ \nabla \psi_h : \psi_h \in W_h \subset H^1 \} \]

Let $u_h \in V_0$, and $\psi_h$ such that $\nabla \psi_h = u_h$. Now, decompose $\psi_h = \sum \psi_i$. Then, $\sum \nabla \psi_i$ is a decomposition of $u_h$ into null-space functions $u_i = \nabla \psi_i$. The block-Jacobi must consist of blocks large enough such that each component $\nabla \psi_i$ is contained in one block. If $\psi_i$ is a vertex-hat basis function, then its gradient is contained in the vertex patch. These block-smoothers are called AFW-blocks, due to Arnold, Falk, and Winther: Multigrid in H(div) and H(curl), Numerische Mathematik, 2000.

Lets try:

from ngsolve import *
from netgen.csg import *
from ngsolve.webgui import Draw

geo = CSGeometry()
magnet = Cylinder(Pnt(-1,0,0), Pnt(1,0,0), 0.3) \
    * OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
box = OrthoBrick(Pnt(-3, -3, -3), Pnt(3,3,3))
geo.Add (magnet.mat("magnet"))
geo.Add ((box-magnet).mat("air"))
mesh = Mesh(geo.GenerateMesh(maxh=1))
# mesh.ngmesh.Refine()

clipping = { "pnt" : (0,0,0), "function" : False}
Draw (mesh, clipping=clipping);

fes = HCurl(mesh, order=0)   # one dof per edge
u,v = fes.TnT()
eps = 1e-3

a = BilinearForm(curl(u)*curl(v)*dx + eps*u*v*dx).Assemble()
M = mesh.MaterialCF( { "magnet" : (1,0,0)} )
f = LinearForm(M*curl(v)*dx("magnet")).Assemble()

gfu = GridFunction(fes)

pre = Preconditioner(a, "local", block=True)
pre.Update()
from ngsolve.krylovspace import CGSolver

inv = CGSolver(a.mat, pre.mat, printrates=True, maxiter=200)
gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.46707005230506027     
CG iteration 2, residual = 0.4201272676774075     
CG iteration 3, residual = 0.31005013098059253     
CG iteration 4, residual = 0.26359732327535396     
CG iteration 5, residual = 0.19697122132018022     
CG iteration 6, residual = 0.15656354546988188     
CG iteration 7, residual = 0.128599226962556     
CG iteration 8, residual = 0.09979601407637312     
CG iteration 9, residual = 0.08090310448574983     
CG iteration 10, residual = 0.06674639074913877     
CG iteration 11, residual = 0.05441920760857872     
CG iteration 12, residual = 0.04495936895306964     
CG iteration 13, residual = 0.038091431470802475     
CG iteration 14, residual = 0.033220434571228294     
CG iteration 15, residual = 0.029735786800423887     
CG iteration 16, residual = 0.02493679288837524     
CG iteration 17, residual = 0.018868517080170726     
CG iteration 18, residual = 0.013705027005993178     
CG iteration 19, residual = 0.009879926829643509     
CG iteration 20, residual = 0.007049761637240523     
CG iteration 21, residual = 0.004839963653361011     
CG iteration 22, residual = 0.0033024654641654145     
CG iteration 23, residual = 0.0022759915680779777     
CG iteration 24, residual = 0.0016085774948354937     
CG iteration 25, residual = 0.0011595466576443806     
CG iteration 26, residual = 0.0008534167553204677     
CG iteration 27, residual = 0.0006994341876678081     
CG iteration 28, residual = 0.0005979781320340812     
CG iteration 29, residual = 0.0004955028723370319     
CG iteration 30, residual = 0.0003969505548331124     
CG iteration 31, residual = 0.00030631216092253446     
CG iteration 32, residual = 0.00022277587161359714     
CG iteration 33, residual = 0.0001656780719595367     
CG iteration 34, residual = 0.000137894492061252     
CG iteration 35, residual = 0.00012153671445432063     
CG iteration 36, residual = 0.00010679095058754028     
CG iteration 37, residual = 9.325598193283489e-05     
CG iteration 38, residual = 7.958067245057104e-05     
CG iteration 39, residual = 6.448452919486306e-05     
CG iteration 40, residual = 5.081946390671325e-05     
CG iteration 41, residual = 4.1426479137999045e-05     
CG iteration 42, residual = 3.432983131123009e-05     
CG iteration 43, residual = 2.902196694522802e-05     
CG iteration 44, residual = 2.501445135273498e-05     
CG iteration 45, residual = 2.0773435784565185e-05     
CG iteration 46, residual = 1.686499157400771e-05     
CG iteration 47, residual = 1.4207337085888709e-05     
CG iteration 48, residual = 1.2585159658073797e-05     
CG iteration 49, residual = 1.169686235456271e-05     
CG iteration 50, residual = 1.0851719643991632e-05     
CG iteration 51, residual = 9.833373837023148e-06     
CG iteration 52, residual = 8.90507314551109e-06     
CG iteration 53, residual = 8.213527009514e-06     
CG iteration 54, residual = 8.001405245978005e-06     
CG iteration 55, residual = 8.020490914782642e-06     
CG iteration 56, residual = 7.852922840440824e-06     
CG iteration 57, residual = 7.380884486271815e-06     
CG iteration 58, residual = 6.8931107805497335e-06     
CG iteration 59, residual = 6.205861847315123e-06     
CG iteration 60, residual = 5.449782861299811e-06     
CG iteration 61, residual = 4.762634313010985e-06     
CG iteration 62, residual = 3.878270476209618e-06     
CG iteration 63, residual = 2.8764184391321674e-06     
CG iteration 64, residual = 2.0047388495914567e-06     
CG iteration 65, residual = 1.4795943738352168e-06     
CG iteration 66, residual = 1.1722527990553722e-06     
CG iteration 67, residual = 1.0034632041040169e-06     
CG iteration 68, residual = 9.124430218716636e-07     
CG iteration 69, residual = 8.380228478309301e-07     
CG iteration 70, residual = 8.088662144708138e-07     
CG iteration 71, residual = 8.155214501533842e-07     
CG iteration 72, residual = 8.097383265449364e-07     
CG iteration 73, residual = 7.802223140271087e-07     
CG iteration 74, residual = 7.494218249536553e-07     
CG iteration 75, residual = 7.177505897095188e-07     
CG iteration 76, residual = 6.535798274655384e-07     
CG iteration 77, residual = 6.131686124648945e-07     
CG iteration 78, residual = 5.608345952057626e-07     
CG iteration 79, residual = 4.791337167005972e-07     
CG iteration 80, residual = 3.8436695533901494e-07     
CG iteration 81, residual = 2.914405720004661e-07     
CG iteration 82, residual = 2.207271017987613e-07     
CG iteration 83, residual = 1.6972547680035683e-07     
CG iteration 84, residual = 1.3693528870687468e-07     
CG iteration 85, residual = 1.1424076605585269e-07     
CG iteration 86, residual = 9.931058462086764e-08     
CG iteration 87, residual = 8.750324133014213e-08     
CG iteration 88, residual = 7.748756487212759e-08     
CG iteration 89, residual = 6.639359879703511e-08     
CG iteration 90, residual = 5.478318267153566e-08     
CG iteration 91, residual = 4.527169952494847e-08     
CG iteration 92, residual = 3.8013208450319486e-08     
CG iteration 93, residual = 3.1441266829476405e-08     
CG iteration 94, residual = 2.5342274064345143e-08     
CG iteration 95, residual = 2.028072308463784e-08     
CG iteration 96, residual = 1.748467100192783e-08     
CG iteration 97, residual = 1.6293418561894208e-08     
CG iteration 98, residual = 1.526996842458048e-08     
CG iteration 99, residual = 1.4335328986167169e-08     
CG iteration 100, residual = 1.2666122069154087e-08     
CG iteration 101, residual = 1.0605004545305719e-08     
CG iteration 102, residual = 8.84921894213243e-09     
CG iteration 103, residual = 7.365771148952591e-09     
CG iteration 104, residual = 6.3265936956471545e-09     
CG iteration 105, residual = 5.708739312109303e-09     
CG iteration 106, residual = 5.213052729288675e-09     
CG iteration 107, residual = 4.7049259181295435e-09     
CG iteration 108, residual = 4.298578835619068e-09     
CG iteration 109, residual = 3.9389300737897795e-09     
CG iteration 110, residual = 3.6975923136629067e-09     
CG iteration 111, residual = 3.485693452586054e-09     
CG iteration 112, residual = 3.2166076956843248e-09     
CG iteration 113, residual = 2.941731890168267e-09     
CG iteration 114, residual = 2.4425853198149036e-09     
CG iteration 115, residual = 1.9483419825102872e-09     
CG iteration 116, residual = 1.525922344964883e-09     
CG iteration 117, residual = 1.186779007656298e-09     
CG iteration 118, residual = 9.445387108893029e-10     
CG iteration 119, residual = 7.805282150747267e-10     
CG iteration 120, residual = 6.207269984797872e-10     
CG iteration 121, residual = 4.820835011451272e-10     
CG iteration 122, residual = 3.906138912797948e-10     
CG iteration 123, residual = 3.4321556643268636e-10     
CG iteration 124, residual = 3.2474317781205145e-10     
CG iteration 125, residual = 3.014558530780509e-10     
CG iteration 126, residual = 2.7255348988950975e-10     
CG iteration 127, residual = 2.3613738098056775e-10     
CG iteration 128, residual = 1.9685930526606346e-10     
CG iteration 129, residual = 1.518978062541263e-10     
CG iteration 130, residual = 1.1294685603637539e-10     
CG iteration 131, residual = 8.50771693975358e-11     
CG iteration 132, residual = 6.582356187810025e-11     
CG iteration 133, residual = 5.416522063962297e-11     
CG iteration 134, residual = 4.619152679951917e-11     
CG iteration 135, residual = 4.008166613611428e-11     
CG iteration 136, residual = 3.6341253046528186e-11     
CG iteration 137, residual = 3.449409329693925e-11     
CG iteration 138, residual = 3.328759415184094e-11     
CG iteration 139, residual = 3.1273244838208485e-11     
CG iteration 140, residual = 2.7747296054761913e-11     
CG iteration 141, residual = 2.4039508476243692e-11     
CG iteration 142, residual = 2.1738026966108458e-11     
CG iteration 143, residual = 2.038848828359549e-11     
CG iteration 144, residual = 1.8496078877749687e-11     
CG iteration 145, residual = 1.6383402924408e-11     
CG iteration 146, residual = 1.4074151489805258e-11     
CG iteration 147, residual = 1.1669667686938074e-11     
CG iteration 148, residual = 9.835363185241042e-12     
CG iteration 149, residual = 8.712281975943535e-12     
CG iteration 150, residual = 8.28327167625019e-12     
CG iteration 151, residual = 7.796147527099701e-12     
CG iteration 152, residual = 7.037248730959671e-12     
CG iteration 153, residual = 6.018744475045147e-12     
CG iteration 154, residual = 5.2519133060422295e-12     
CG iteration 155, residual = 4.784599186787955e-12     
CG iteration 156, residual = 4.200866793083326e-12     
CG iteration 157, residual = 3.5544804597975956e-12     
CG iteration 158, residual = 3.0079033517655157e-12     
CG iteration 159, residual = 2.52325188281144e-12     
CG iteration 160, residual = 2.0899754351058183e-12     
CG iteration 161, residual = 1.6912229341695938e-12     
CG iteration 162, residual = 1.4045133157463545e-12     
CG iteration 163, residual = 1.1720721649793027e-12     
CG iteration 164, residual = 9.647111470300315e-13     
CG iteration 165, residual = 7.824611455164816e-13     
CG iteration 166, residual = 6.233659600891288e-13     
CG iteration 167, residual = 5.228207325064967e-13     
CG iteration 168, residual = 4.433011957647366e-13     

Draw (curl(gfu), mesh, clipping=clipping, \
      vectors={"grid_size" : 40 }, draw_surf=False);

from ngsolve.la import EigenValues_Preconditioner
lam = EigenValues_Preconditioner(a.mat, pre.mat)
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 0.02939700558650116 6.141387994025265
kappa= 208.91202595292

We observe that the condition number of the block-Jacobi preconditioner

is robust in $\varepsilon$
depends on the mesh-size

We observe the the condition number of the Jacobi preconditioner

depends on $\varepsilon$
depends on the mesh-size

57.2. Robust coarse-grid correction#

If the penalty term involves a projection, then the coarse-grid null-space $V_{H,0}$ is not a sub-space of the fine-grid null-space $V_{h,0}$. An example is Stokes with penalty:

\[ \| u_H \|_{A_H^\varepsilon}^2 = \| u_H \|_A^2 + \frac{1}{\varepsilon} \| R_H \operatorname{div} u_H \|_{L_2}^2 \]

The operator-norm of the prolongation is $$ \| P \|_{(V_H, A_H^\varepsilon) \rightarrow (V_h, A_h^\varepsilon)} = \sup_{v_H} \frac{ \| P v_H \|_{A_h^\varepsilon}}{\| v_H \|_{A_H^\varepsilon}} $$

To be robust in the parameter, the prolongation operator $P : V_H \rightarrow V_h$ must map the coarse-grid null-space into the fine-grid null-space:

\[ P V_{H,0} \subset V_{h,0} \]

If the null-spaces are nested (as for Maxwell equations), the prolongation is canonical embedding.

For the ASM - analysis we need the existence of stable interpolation operators from fine to coarse: $I_H : V_h \rightarrow V_H$. It must be compatible with null-spaces, i.e. $I_H V_{h,0} \subset V_{H,0}$.

57.3. Two-level analysis for Maxwell equations#

Let $V_h$ be a Nedelec finite element sub-space of $H(\operatorname{curl})$, and consider the bilinear-form

\[ A^\varepsilon(u,v) = \int \operatorname{curl} u \operatorname{curl} v + \varepsilon u v \]

We prove that a two-level method with an AFW block-Jacobi smoother provides an optimal preconditioner.

Let $W_H \subset W_h \subset H^1$ be compatible finite element spaces such that $\nabla W_h = V_{h,0} = \{ v \in V_h : \operatorname{curl} v_h = 0 \}$, and $\nabla W_H = W_{H,0}$.

We need a few technical tools:

Theorem: (Regular Decomposition): There holds

\[ H(\operatorname{curl}) = [H^1]^3 + \nabla H^1 \]

For every function $u$ from $H(\operatorname{curl})$ there exists a decomposition

\[ u = z + \nabla \phi \]

with $z \in [H^1]^3$ and $\phi \in H^1$, and which is stable in the sense

\[\begin{eqnarray*} \| \phi \|_{H^1} & \preceq & \| u \|_{H(\operatorname{curl})} \\ \| z \|_{H^1} & \preceq & \| \operatorname{curl} u \|_{L_2} \\ \end{eqnarray*}\]

On the other hand, every function from $[H^1]^3$, and every gradient from $H^1$ is in $H(\operatorname{curl})$.

Proof: Pasciak, J. E.; Zhao, J.: Overlapping Schwarz methods in H(curl) on polyhedral domains. J. Numer. Math. 10 (2002)

Theorem: (Commuting quasi-interpolation operators) There exist projections $\pi_{V_h} : H(\operatorname{curl}) \rightarrow V_h $ and $\pi_{W_h} : H^1 \rightarrow W_h$ which commute with the gradient:

\[ \pi_{V_h} \nabla = \nabla \pi_{W_h} \]

They are stable in norms and semi-norms, and provide optimal approximation error estimates.

Proof: Schöberl, A multilevel decomposition result in H(curl). in Multigrid, Multilevel and Multiscale Methods, EMG 2005

We are ready to prove that the condition-number of the two-level preconditioner is robust in the mesh-size $h$, as well as the regularization parameter $\varepsilon$. We assume that $H = ch$. We transfer results from the $H^1$-case to the $H(\operatorname{curl})$ case.

In the $H^1$-case, we have shown that the coarse grid interpolation error satisfies

\[ \| u_h - \pi_{W_H} u_h \|_{L_2} \preceq h \| \nabla u_h \| \]

The corresponding lemma for the $H(\operatorname{curl})$ case is:

Lemma: For the coarse-grid interpolation error there exists a decomposition

\[ u_h - \pi_{V_H} u_h = r + \nabla \psi, \]

where $$ \| r \|_{L_2}^2 \preceq h^2 \| \operatorname{curl} u \|^2 $$ and $$ \varepsilon \| \psi \|_{L_2}^2 \preceq h^2 \| u \|_{A^\varepsilon}^2 $$

Proof: Let $u_h \in V_h \subset H(\operatorname{curl})$. There exists a regular decomposition

\[ u_h = z + \nabla \Phi, \]

where

\[ \| z \|_{H^1}^2 \preceq \| \operatorname{curl} u \|_{L_2}^2 \]

and $$ \varepsilon \| \Phi \|_{H^1}^2 \preceq \| u \|_{A^\varepsilon}^2 $$

By the commutativity of the interpolation operator we obtain

\[\begin{eqnarray*} u_h - \pi_{V_H} u_h & = & z - \pi_{V_H} z + \nabla \Phi - \pi_{V_H} \nabla \Phi \\ & = & z - \pi_{V_H} z + \nabla ( \Phi - \pi_{W_H} \Phi ) \end{eqnarray*}\]

By setting

\[ r = z - \pi_{V_H} z \qquad \text{and} \qquad \psi = \Phi - \pi_{W_H} \Phi \]

and using the approximation properties of the interpolation operators we obtain the splitting.

Lemma: For the coarse-grid interpolation error there exists a discrete decomposition

\[ u_h - \pi_{V_H} u_h = r_h + \nabla \psi_h. \]

Proof: Apply the commuting projection operator $\pi_{V_h}$ to the previous lemma to obtain

(57.1)#\[\begin{eqnarray} u_h - \pi_{V_H} u_h & = & \pi_{V_h} ( u_h - \pi_{V_H} u_h ) \\ & = & \pi_{V_h} r + \pi_{V_h} \nabla \psi \\ & = & \pi_{V_h} r + \nabla \pi_{W_h} \psi \end{eqnarray}\]

The discrete decomposition is now $r_h = \pi_{V_h} r$ and $\psi_h = \pi_{W_h} \psi$.

In the $H^1$-case, we can apply the ASM lemma to analyze the Jacobi preconditioner $D$ as follows:

\[ \| u \|_D^2 = \inf_{u = \sum u_i} \sum \| u_i \|_{H^1}^2 \preceq \inf_{u = \sum u_i} h^{-2} \sum \| u_i \|_{L_2}^2 \preceq h^{-2} \| u \|_{L_2}^2 \]

Lemma: For the AFW block-Jacobi preconditioner $D$ there holds

\[ \| u_h \|_D^2 \preceq \inf_{u_h = r_h + \nabla \psi_h} h^{-2} \| r_h \|_{L_2}^2 + \varepsilon h^{-2} \| \psi_h \|_{L_2}^2 \]

Proof: Let $u_h = r_h + \nabla \psi_h$ be an arbitrary decomposition. Then the triangle inequality implies

\[ \| u_h \|_D \leq \| r_h \|_D + \| \nabla \psi_h \|_D \]

Now we apply the ASM lemma for both terms:

\[\begin{eqnarray*} \| r_h \|_D^2 & = & \inf_{r = \sum r_i} \sum \| r_i \|_{A^\varepsilon}^2 \\ & = & \inf_{r = \sum r_i} \sum \| \operatorname{curl} r_i \|_{L_2}^2 + \varepsilon \| r_i \|_{L_2}^2 \\ & \preceq & \inf_{r = \sum r_i} h^{-2} \| r_i \|_{L_2}^2 \\ & \approx & h^{-2} \| r \|_{L_2}^2 \end{eqnarray*}\]

For the decomposition of $\nabla \psi$ we use that we can decompose $\psi = \sum \psi_i$ with $\psi_i \in W_i$, where $\nabla W_i \subset V_i$:

\[\begin{eqnarray*} \varepsilon \| \nabla \psi_h \|_D^2 & = & \inf_{\nabla \psi_h = \sum v_i} \sum \| v_i \|_{A^\varepsilon}^2 \\ & \leq & \inf_{\psi_h = \sum \psi_i} \sum \| \nabla \psi_i \|_{A^\varepsilon}^2 \\ & = & \inf_{\psi_h = \sum \psi_i} \sum \varepsilon \| \nabla \psi_i \|_{L_2}^2 \\ & \preceq & \inf_{\psi_h = \sum \psi_i} \sum \varepsilon h^{-2} \| \psi_i \|_{L_2}^2 \\ & \approx & \varepsilon h^{-2} \| \psi_h \|_{L_2}^2 \end{eqnarray*}\]

Since $\| u \|_D$ is bounded by any such decomposition of $u_h$, it holds also for the infimum, and the lemma is proven.

Theorem The two-level preconditioner with AFW block-smoothers provide a optimal condition number $\kappa = O(1)$.

Proof: We have to verify that

\[ \inf_{u_h = u_H + \sum u_i} \left\{ \| u_H \|_{A^\varepsilon}^2 + \sum \| u_i \|_{A^\varepsilon}^2 \right\} \preceq \| u_h \|_{A^\varepsilon}^2 \]

We choose $u_H := \pi_{V_H} u_h$. Continuity in energy-norm of the projector bounds the first term. The interpolation rest provides a stable splitting

\[ u_f = u_h - u_H = r_h + \nabla \psi_h, \]

where the components $r_h$ and $\psi_h$ are bounded in $L_2$-norms. This $u_f$ can be further split into local functions as $u_f = \sum u_i$.

from ngsolve import *
from netgen.csg import *
from ngsolve.webgui import Draw

geo = CSGeometry()
magnet = Cylinder(Pnt(-1,0,0), Pnt(1,0,0), 0.3) \
    * OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
box = OrthoBrick(Pnt(-3, -3, -3), Pnt(3,3,3))
geo.Add (magnet.mat("magnet"))
geo.Add ((box-magnet).mat("air"))
mesh = Mesh(geo.GenerateMesh(maxh=1))

fes = HCurl(mesh, order=0, autoupdate=True)
print ("ndof =", fes.ndof)

u,v = fes.TnT()
eps = 1e-3

a = BilinearForm(curl(u)*curl(v)*dx + eps*u*v*dx)
pre = Preconditioner(a, "multigrid", smoother="block")

M = mesh.MaterialCF( { "magnet" : (1,0,0)} )
f = LinearForm(M*curl(v)*dx("magnet"))

with TaskManager():
    a.Assemble()
    f.Assemble()
    for l in range(1):  # try also 2
        mesh.Refine()
        print ("ndof =", fes.ndof)
        a.Assemble()
        f.Assemble()

gfu = GridFunction(fes)

ndof = 10124

WARNING: kwarg 'smoother' is an undocumented flags option for class <class 'ngsolve.comp.Preconditioner'>, maybe there is a typo?

ndof = 89289

from ngsolve.krylovspace import CGSolver

inv = CGSolver(a.mat, pre.mat, printrates=True, maxiter=200)
with TaskManager():
    gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.6907924751199783     
CG iteration 2, residual = 0.01695481662307981     
CG iteration 3, residual = 0.0005294588214023937     

CG iteration 4, residual = 1.2097981420766042e-05     
CG iteration 5, residual = 3.860394989844051e-07     
CG iteration 6, residual = 3.093473387425276e-08     

CG iteration 7, residual = 2.5824941848269288e-09     
CG iteration 8, residual = 1.3354858061212555e-10     
CG iteration 9, residual = 1.3782726321179078e-11     

CG iteration 10, residual = 6.847534975421774e-13

clipping = { "pnt" : (0,0,0), "function" : False}
Draw (curl(gfu), mesh, clipping=clipping, \
      vectors={"grid_size" : 40 }, draw_surf=False);